{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Surface restructuring event analysis of large-scale LAMMPS MD trajectory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Header"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import lmpdump as lmpdump    # Use lmpdump.py parser\n",
    "import os\n",
    "import math\n",
    "import numpy as np\n",
    "import time\n",
    "from tqdm import tqdm\n",
    "\n",
    "pwd = os.getcwd()\n",
    "\n",
    "file_out = pwd + '/output.txt'\n",
    "sys.stdout = open(file_out, 'w')\n",
    "\n",
    "print('Please refer to README for more information.\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read input file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_in = pwd + '/input.txt'\n",
    "Input = open(file_in,'r')\n",
    "input_lines = Input.readlines()\n",
    "input_lines = [line.split() for line in input_lines]\n",
    "\n",
    "for line in input_lines:\n",
    "    if line[0] != '#':\n",
    "    \n",
    "        if 'clamp_dump' in line:\n",
    "            clamp_dump = line[1]\n",
    "\n",
    "        elif 'eq_data' in line:\n",
    "            eq_data = line[1]\n",
    "\n",
    "        elif 'ID_NN1' in line:\n",
    "            ID_nn1 = int(line[1])\n",
    "\n",
    "        elif 'ID_NN2' in line:\n",
    "            ID_nn2 = int(line[1])\n",
    "        \n",
    "        elif 'ID_1up' in line:\n",
    "            ID_1up = int(line[1])\n",
    "\n",
    "if (clamp_dump and eq_data and ID_nn1 and ID_nn2 and ID_1up):\n",
    "    print('Done reading input file.\\n')\n",
    "else:\n",
    "    print('Error: Input variable missing.\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load LAMMPS trajectory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# LAMMPS custom-style dump file - Clamped trajectory\n",
    "print('Loading trajectory...')\n",
    "file = pwd + '/' + clamp_dump\n",
    "xyz = lmpdump.lmpdump(file, loadmode='all')\n",
    "print('\\n')\n",
    "\n",
    "print('Loading time steps...')\n",
    "step = []\n",
    "for s in xyz.finaldict.keys():    # keys = time steps\n",
    "    step.append(s)\n",
    "print('Done.\\n')\n",
    "\n",
    "size = np.array(step).size    # Number of time steps\n",
    "N_dump = xyz.finaldict[0][1].shape[0]    # Number of atoms dumped\n",
    "xlo_ext = xyz.finaldict[0][0][0]    # x_ext = x + xy\n",
    "xhi_ext = xyz.finaldict[0][0][1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract unit cell information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load LAMMPS data file after NPT equilibration\n",
    "file_eq = pwd + '/' + eq_data\n",
    "log = open(file_eq,'r')\n",
    "log_lines = log.readlines()\n",
    "log_lines = [line.split() for line in log_lines]\n",
    "\n",
    "for line_index, line in enumerate(log_lines):\n",
    "    if line:\n",
    "        for l in line:\n",
    "            \n",
    "            if l == 'atoms':\n",
    "                N_atom = int(line[0])\n",
    "                \n",
    "            elif l == 'types':\n",
    "                N_type = int(line[0])\n",
    "\n",
    "            elif l == 'xlo':\n",
    "                xlo = float(line[0])\n",
    "                xhi = float(line[1])\n",
    "            \n",
    "            elif l == 'ylo':\n",
    "                ylo = float(line[0])\n",
    "                yhi = float(line[1])\n",
    "            \n",
    "            elif l == 'zlo':\n",
    "                zlo = float(line[0])\n",
    "                zhi = float(line[1])\n",
    "            \n",
    "            elif l == 'xy':\n",
    "                xy = float(line[0])\n",
    "                xz = float(line[1])\n",
    "                yz = float(line[2])\n",
    "            \n",
    "            elif l == 'Atoms':\n",
    "                head = line_index + 1\n",
    "                start = line_index + 2\n",
    "            \n",
    "            elif l == 'Velocities':\n",
    "                end = line_index - 1\n",
    "\n",
    "Lx = xhi - xlo    # Box dimensions\n",
    "Ly = yhi - ylo\n",
    "Lz = zhi - zlo\n",
    "tan = Ly/xy\n",
    "\n",
    "xyz = []    # N*3 matrix\n",
    "for line in log_lines[start:end]:\n",
    "\n",
    "    line[:2] = np.array(line[:2]).astype(int)    # Atom ID, atom type\n",
    "    line[2:5] = np.array(line[2:5]).astype(float)    # x, y, z\n",
    "    line[5:8] = np.array(line[5:8]).astype(int)    # ix, iy, iz (inverse sign)\n",
    "\n",
    "    x = line[2] + line[5]*Lx + line[6]*xy    # Wrapped coordinates\n",
    "    y = line[3] + line[6]*Ly\n",
    "    z = line[4] + line[7]*Lz\n",
    "    \n",
    "    xyz.append([x,y,z])\n",
    "xyz = np.array(xyz)\n",
    "\n",
    "ID2 = int(input('Enter the ID of a nearest-neighbor atom that lies in the same layer as atom #1 : '))\n",
    "ID3 = int(input('Enter the ID of another nearest-neighbor atom that lies in the same layer as atom #1 : '))\n",
    "ID4 = int(input('Enter the ID of an atom that lies one layer above atom #1 : '))\n",
    "\n",
    "for atom, coord in enumerate(xyz):\n",
    "    \n",
    "    if atom+1 == 1:\n",
    "        atom1 = coord    # atom1 = Reference atom\n",
    "    \n",
    "    elif atom+1 == ID2:\n",
    "        atom2 = coord    # atom2 = 1st coplanar atom\n",
    "\n",
    "    elif atom+1 == ID3:\n",
    "        atom3 = coord    # atom3 = 2nd coplanar atom\n",
    "        \n",
    "    elif atom+1 == ID4:\n",
    "        atom4 = coord    # atom4 = Atom one layer above\n",
    "\n",
    "r12 = np.linalg.norm(np.subtract(atom1,atom2))\n",
    "r13 = np.linalg.norm(np.subtract(atom1,atom3))\n",
    "dr_avg = np.mean([r12, r13])    # In-plane NN distance\n",
    "\n",
    "a1 = np.subtract(atom2,atom1)    # 1st intralayer vector\n",
    "a1 = a1 / np.linalg.norm(a1)\n",
    "\n",
    "a2 = np.subtract(atom3,atom1)    # 2nd intralayer vector\n",
    "a2 = a2 / np.linalg.norm(a2)\n",
    "\n",
    "a3 = np.cross(a1,a2)    # Normal vector\n",
    "if a3[2] < 0:\n",
    "    a3 = -a3    # Ensure normal vector along +z\n",
    "a3 = a3 / np.linalg.norm(a3)\n",
    "\n",
    "dn_avg = np.dot(atom4,a3) - np.dot(atom1,a3)    # Interlayer distance; n = normal component\n",
    "\n",
    "N_slab = int(input('Enter the number of layers in the slab, excluding the deposit/adsorbate layer(s) : '))\n",
    "if xyz1.finaldict[0][1]['id'][0] != 1:    # Is bottommost layer muted?\n",
    "    N_slab -= 1\n",
    "\n",
    "z_layer = []    # Clamped layer heights\n",
    "for atom in range(1,N_dump1):\n",
    "    z = xyz1.finaldict[step[-1]][1]['z'][atom]    # Use last frame\n",
    "    if z not in z_layer:\n",
    "        z_layer.append(z)\n",
    "N_layer = len(z_layer)\n",
    "\n",
    "log.close()\n",
    "\n",
    "print('Unit cell information & normal direction extracted.\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Event detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Parsing trajectory for restructuring events...')\n",
    "t_begin = time.perf_counter()\n",
    "\n",
    "file_ev = pwd + '/event_flag.txt'    # List of time steps flagged as event points\n",
    "out_ev = open(file_ev,'w')\n",
    "\n",
    "header1 = ' tf = final time (flagged step) \\n dt = event duration \\n ni = initial normal component \\n nf = final normal \\n CNi = initial coordination number \\n CNf = final coordination number \\n \\n'\n",
    "header2 = 'tf (ps) \\t dt (ps) \\t ID \\t type \\t ni \\t nf \\t CNi \\t CNf \\n'\n",
    "out_ev.write(header1)\n",
    "out_ev.write(header2)\n",
    "\n",
    "skip_index = [0, 1, 2, size-2, size-1]    # Skip time steps at the edges of the trajectory\n",
    "\n",
    "for s_index, s in enumerate(tqdm(step)):\n",
    "    \n",
    "    if s_index in skip_index:\n",
    "        continue\n",
    "\n",
    "    ev_atom = []    # List of active atoms detected\n",
    "    for atom in range(0,N_dump):\n",
    "            \n",
    "        # f0 = current time step (final 0)\n",
    "        # f1 = 1 step after (final 1)\n",
    "        # f2 = 2 steps after (final 2)\n",
    "        # i1 = 1 step before (initial 1)\n",
    "        # i2 = 2 steps before (initial 2)\n",
    "        # i3 = 3 steps before (initial 3)\n",
    "        \n",
    "        ID = xyz1.finaldict[s][1]['id'][atom]\n",
    "        Type = xyz1.finaldict[s][1]['type'][atom]\n",
    "        \n",
    "        CNf0 = xyz.finaldict[s][1]['c_cn'][atom]\n",
    "        CNi1 = xyz.finaldict[step[s_index-1]][1]['c_cn'][atom]\n",
    "        CNi2 = xyz.finaldict[step[s_index-2]][1]['c_cn'][atom]\n",
    "        \n",
    "        xf0 = xyz.finaldict[s][1]['xu'][atom]\n",
    "        yf0 = xyz.finaldict[s][1]['yu'][atom]\n",
    "        zf0 = xyz.finaldict[s][1]['z'][atom]\n",
    "        rf0 = np.array([xf0, yf0, zf0])\n",
    "        nf0 = np.dot(rf0,a3)\n",
    "        xf0_prj = xf0 - yf0/tan    # Tilted projection of x\n",
    "        Nxf0 = math.floor(xf0_prj/Lx)    # x image number\n",
    "        Nyf0 = math.floor(yf0/Ly)    # y image number\n",
    "    \n",
    "        xf1 = xyz.finaldict[step[s_index+1]][1]['xu'][atom]\n",
    "        yf1 = xyz.finaldict[step[s_index+1]][1]['yu'][atom]\n",
    "        zf1 = xyz.finaldict[step[s_index+1]][1]['z'][atom]\n",
    "        rf1 = np.array([xf1, yf1, zf1])\n",
    "        nf1 = np.dot(rf1,a3)\n",
    "        xf1_prj = xf1 - yf1/tan\n",
    "        Nxf1 = math.floor(xf1_prj/Lx)\n",
    "        Nyf1 = math.floor(yf1/Ly)\n",
    "        \n",
    "        xf2 = xyz.finaldict[step[s_index+2]][1]['xu'][atom]\n",
    "        yf2 = xyz.finaldict[step[s_index+2]][1]['yu'][atom]\n",
    "        zf2 = xyz.finaldict[step[s_index+2]][1]['z'][atom]\n",
    "        rf2 = np.array([xf2, yf2, zf2])\n",
    "        nf2 = np.dot(rf2,a3)\n",
    "        xf2_prj = xf2 - yf2/tan\n",
    "        Nxf2 = math.floor(xf2_prj/Lx)\n",
    "        Nyf2 = math.floor(yf2/Ly)\n",
    "        \n",
    "        xi1 = xyz.finaldict[step[s_index-1]][1]['xu'][atom]\n",
    "        yi1 = xyz.finaldict[step[s_index-1]][1]['yu'][atom]\n",
    "        zi1 = xyz.finaldict[step[s_index-1]][1]['z'][atom]\n",
    "        ri1 = np.array([xi1, yi1, zi1])\n",
    "        ni1 = np.dot(ri1,a3)\n",
    "        xi1_prj = xi1 - yi1/tan\n",
    "        Nxi1 = math.floor(xi1_prj/Lx)\n",
    "        Nyi1 = math.floor(yi1/Ly)\n",
    "        \n",
    "        xi2 = xyz.finaldict[step[s_index-2]][1]['xu'][atom]\n",
    "        yi2 = xyz.finaldict[step[s_index-2]][1]['yu'][atom]\n",
    "        zi2 = xyz.finaldict[step[s_index-2]][1]['z'][atom]\n",
    "        ri2 = np.array([xi2, yi2, zi2])\n",
    "        ni2 = np.dot(ri2,a3)\n",
    "        xi2_prj = xi2 - yi2/tan\n",
    "        Nxi2 = math.floor(xi2_prj/Lx)\n",
    "        Nyi2 = math.floor(yi2/Ly)\n",
    "        \n",
    "        xi3 = xyz.finaldict[step[s_index-3]][1]['xu'][atom]\n",
    "        yi3 = xyz.finaldict[step[s_index-3]][1]['yu'][atom]\n",
    "        zi3 = xyz.finaldict[step[s_index-3]][1]['z'][atom]\n",
    "        ri3 = np.array([xi3, yi3, zi3])\n",
    "        ni3 = np.dot(ri3,a3)\n",
    "        xi3_prj = xi3 - yi3/tan\n",
    "        Nxi3 = math.floor(xi3_prj/Lx)\n",
    "        Nyi3 = math.floor(yi3/Ly)\n",
    "    \n",
    "        dni1 = nf0 - ni1    # Change in normal from the past\n",
    "        dni2 = nf0 - ni2\n",
    "        dni3 = nf0 - ni3\n",
    "        \n",
    "        dnf1 = nf1 - nf0    # Change in normal in the future\n",
    "        dnf2 = nf2 - nf0\n",
    "        \n",
    "        # Test 1 step in the past\n",
    "        # Criterion: Change by more than 90% of interlayer distance\n",
    "        # History check: Must not revert back 1 or 2 steps in the past & in the future\n",
    "        if (abs(dni1) > 0.90*dn_avg) and (abs(dni2) > 0.90*dn_avg) and (abs(dnf1) < 0.10*dn_avg) and (abs(dnf2) < 0.10*dn_avg):\n",
    "            \n",
    "            ev_atom.append(atom)\n",
    "            dt = 1\n",
    "            \n",
    "            xf0 -= xy * Nyf0    # Wrap y first\n",
    "            yf0 -= Ly * Nyf0\n",
    "            \n",
    "            xi1 -= xy * Nyi1\n",
    "            yi1 -= Ly * Nyi1\n",
    "            \n",
    "            if (Nxf0 != 0) or (Nxi1 != 0):\n",
    "            \n",
    "                xf0 -= Lx * np.maximum(Nxf0,Nxi1)    # Shift initial & final x together toward the unit cell\n",
    "                rf0 = np.array([xf0, yf0, zf0])\n",
    "                nf0 = np.dot(rf0,a3)\n",
    "            \n",
    "                xi1 -= Lx * np.maximum(Nxf0,Nxi1)\n",
    "                ri1 = np.array([xi1, yi1, zi1])\n",
    "                ni1 = np.dot(ri1,a3)\n",
    "            \n",
    "            line = '%i\\t%i\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(s_index, dt, ID, Type, ni1, nf0, CNi1, CNf0)\n",
    "            out_ev.write(line)\n",
    "            \n",
    "            for aux in range(0,N_dump):    # Detect any auxiliary atom that may move together\n",
    "                if aux not in ev_atom:    # Test only if not flagged before\n",
    "                    \n",
    "                    # Auxiliary atom labeled with \"i\" & \"f\" only\n",
    "                    \n",
    "                    ID = xyz.finaldict[s][1]['id'][aux]\n",
    "                    Type = xyz.finaldict[s][1]['type'][aux]\n",
    "                    \n",
    "                    CNf = xyz.finaldict[s][1]['c_cn'][aux]\n",
    "                    CNi = xyz.finaldict[step[s_index-1]][1]['c_cn'][aux]\n",
    "    \n",
    "                    xf = xyz.finaldict[s][1]['xu'][aux]\n",
    "                    yf = xyz.finaldict[s][1]['yu'][aux]\n",
    "                    zf = xyz.finaldict[s][1]['z'][aux]\n",
    "                    rf = np.array([xf, yf, zf])\n",
    "                    nf = np.dot(rf,a3)\n",
    "                    xf_prj = xf - yf/tan\n",
    "                    Nxf = math.floor(xf_prj/Lx)\n",
    "                    Nyf = math.floor(yf/Ly)\n",
    "                    \n",
    "                    xi = xyz.finaldict[step[s_index-1]][1]['xu'][aux]\n",
    "                    yi = xyz.finaldict[step[s_index-1]][1]['yu'][aux]\n",
    "                    zi = xyz.finaldict[step[s_index-1]][1]['z'][aux]\n",
    "                    ri = np.array([xi, yi, zi])\n",
    "                    ni = np.dot(ri,a3)\n",
    "                    xi_prj = xi - yi/tan\n",
    "                    Nxi = math.floor(xi_prj/Lx)\n",
    "                    Nyi = math.floor(yi/Ly)\n",
    "    \n",
    "                    dx = xf - xi\n",
    "                    dy = yf - yi\n",
    "                    dz = zf - zi                        \n",
    "                    dr = math.sqrt(dx**2 + dy**2 + dz**2)    # Change in position, rather than normal\n",
    "                    \n",
    "                    if dr > 0.80*dr_avg:    # Criterion: Change by more than 80% of the NN distance\n",
    "                        ev_atom.append(aux)\n",
    "                        \n",
    "                        # Auxiliary atom must be wrapped toward the active atom\n",
    "\n",
    "                        xf -= xy * Nyf    # Wrap y first\n",
    "                        yf -= Ly * Nyf\n",
    "                        \n",
    "                        xi -= xy * Nyi\n",
    "                        yi -= Ly * Nyi\n",
    "                        \n",
    "                        dxf = (xf - xf0)/Lx    # Extent of x deviation from the active atom\n",
    "                        dxi = (xi - xi1)/Lx\n",
    "                        \n",
    "                        if (abs(dxf) > 0.70) or (abs(dxi) > 0.70):    # Wrap if x deviation is more than 70% of the unit cell x dimension\n",
    "                                \n",
    "                            xf -= Lx * round(dxf)\n",
    "                            rf = np.array([xf, yf, zf])\n",
    "                            nf = np.dot(rf,a3)\n",
    "                            \n",
    "                            xi -= Lx * round(dxi)\n",
    "                            ri = np.array([xi, yi, zi])\n",
    "                            ni = np.dot(ri,a3)\n",
    "\n",
    "                        line = '%i\\t%i\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(s_index, dt, ID, Type, ni, nf, CNi, CNf)\n",
    "                        out_ev.write(line)\n",
    "        \n",
    "        # If no change 1 step in the past, test 2 steps in the past (change may be gradual)\n",
    "        elif (abs(dni2) > 0.90*dn_avg) and (abs(dni3) > 0.90*dn_avg) and (abs(dnf1) < 0.10*dn_avg) and (abs(dnf2) < 0.10*dn_avg): \n",
    "            \n",
    "            ev_atom.append(atom)\n",
    "            dt = 2\n",
    "            \n",
    "            xf0 -= xy * Nyf0\n",
    "            yf0 -= Ly * Nyf0\n",
    "            \n",
    "            xi2 -= xy * Nyi2\n",
    "            yi2 -= Ly * Nyi2\n",
    "    \n",
    "            if (Nxf0 != 0) or (Nxi2 != 0):\n",
    "            \n",
    "                xf0 -= Lx * np.maximum(Nxf0,Nxi2)\n",
    "                rf0 = np.array([xf0, yf0, zf0])\n",
    "                nf0 = np.dot(rf0,a3)\n",
    "            \n",
    "                xi2 -= Lx * np.maximum(Nxf0,Nxi2)\n",
    "                ri2 = np.array([xi2, yi2, zi2])\n",
    "                ni2 = np.dot(ri2,a3)\n",
    "            \n",
    "            line = '%i\\t%i\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(s_index, dt, ID, Type, ni2, nf0, CNi2, CNf0)\n",
    "            out_ev.write(line)\n",
    "    \n",
    "            for aux in range(0,N_dump):\n",
    "                if aux not in ev_atom:\n",
    "                    \n",
    "                    ID = xyz.finaldict[s][1]['id'][aux]\n",
    "                    Type = xyz.finaldict[s][1]['type'][aux]\n",
    "                    \n",
    "                    CNf = xyz.finaldict[s][1]['c_cn'][aux]\n",
    "                    CNi = xyz.finaldict[step[s_index-1]][1]['c_cn'][aux]\n",
    "    \n",
    "                    xf = xyz.finaldict[s][1]['xu'][aux]                            \n",
    "                    yf = xyz.finaldict[s][1]['yu'][aux]\n",
    "                    zf = xyz.finaldict[s][1]['z'][aux]\n",
    "                    rf = np.array([xf, yf, zf])\n",
    "                    nf = np.dot(rf,a3)\n",
    "                    xf_prj = xf - yf/tan\n",
    "                    Nxf = math.floor(xf_prj/Lx)\n",
    "                    Nyf = math.floor(yf/Ly)\n",
    "      \n",
    "                    xi = xyz.finaldict[step[s_index-2]][1]['xu'][aux]\n",
    "                    yi = xyz.finaldict[step[s_index-2]][1]['yu'][aux]\n",
    "                    zi = xyz.finaldict[step[s_index-2]][1]['z'][aux]\n",
    "                    ri = np.array([xi, yi, zi])\n",
    "                    ni = np.dot(ri,a3)\n",
    "                    xi_prj = xi - yi/tan\n",
    "                    Nxi = math.floor(xi_prj/Lx)\n",
    "                    Nyi = math.floor(yi/Ly)\n",
    "    \n",
    "                    dx = xf - xi\n",
    "                    dy = yf - yi\n",
    "                    dz = zf - zi                        \n",
    "                    dr = math.sqrt(dx**2 + dy**2 + dz**2)\n",
    "                    \n",
    "                    if dr > 0.80*dr_avg:\n",
    "                        ev_atom.append(aux)\n",
    "                        \n",
    "                        xf -= xy * Nyf\n",
    "                        yf -= Ly * Nyf\n",
    "                        \n",
    "                        xi -= xy * Nyi\n",
    "                        yi -= Ly * Nyi\n",
    "                        \n",
    "                        dxf = (xf - xf0)/Lx\n",
    "                        dxi = (xi - xi2)/Lx\n",
    "                        \n",
    "                        if (abs(dxf) > 0.70) or (abs(dxi) > 0.70):\n",
    "\n",
    "                            xf -= Lx * round(dxf)\n",
    "                            rf = np.array([xf, yf, zf])\n",
    "                            nf = np.dot(rf,a3)\n",
    "\n",
    "                            xi -= Lx * round(dxi)\n",
    "                            ri = np.array([xi, yi, zi])\n",
    "                            ni = np.dot(ri,a3)\n",
    "                            \n",
    "                        line = '%i\\t%i\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(s_index, dt, ID, Type, ni, nf, CNi, CNf)\n",
    "                        out_ev.write(line)\n",
    "                            \n",
    "out_ev.close()\n",
    "\n",
    "t_end = time.perf_counter()\n",
    "print('Event detection done in %.2f' % (t_end - t_begin) + 's > event_flag.txt\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Event clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ev = open(file_ev,'r')    # Load list of time steps flagged as event points\n",
    "ev_lines = ev.readlines()\n",
    "ev_lines = [line.split() for line in ev_lines]\n",
    "\n",
    "for line_index, line in enumerate(ev_lines):\n",
    "    if line_index > 7:\n",
    "        line[0:4] = np.array(line[0:4]).astype(int)    # Time step, duration, atom ID, atom type\n",
    "        line[4:6] = np.array(line[4:6]).astype(float)    # Initial normal, final normal\n",
    "        line[6:8] = np.array(line[6:8]).astype(int)    # Initial coordination, final coordination\n",
    "\n",
    "ev.close()\n",
    "\n",
    "file_cluster = pwd + '/event_cluster.txt'    # List of events\n",
    "out_cluster = open(file_cluster,'w')\n",
    "\n",
    "# List active atoms first, followed by event #, final time, and duration\n",
    "header1 = ' tf = final time (flagged step) \\n dt = event duration \\n ni = initial normal component \\n nf = final normal component \\n CNi = initial coordination number \\n CNf = final coordination number \\n \\n'\n",
    "header2 = '\\t ID \\t type \\t ni \\t nf \\t CNi \\t CNf \\n event # \\t tf (ps) \\t dt (ps) \\n \\n'\n",
    "out_cluster.write(header1)\n",
    "out_cluster.write(header2)\n",
    "\n",
    "ev_no = 0\n",
    "ev_atom = []    # List of active atoms for each event\n",
    "ev_flag = []    # List of flags for each event\n",
    "\n",
    "for line_index, line in enumerate(ev_lines):\n",
    "\n",
    "    if line_index < 8:\n",
    "        continue\n",
    "    \n",
    "    elif line_index == 8:\n",
    "        \n",
    "        ev_flag.append(line)\n",
    "        ti = line[0] - line[1]    # ti = initial time\n",
    "        ID = line[2]\n",
    "        ev_atom.append(ID)\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        # If a step is within 3 steps of the previously flagged step, it is part of the same event\n",
    "        if line[0] - ev_lines[line_index-1][0] < 3:\n",
    "            \n",
    "            ev_flag.append(line)\n",
    "            ID = line[2]\n",
    "            \n",
    "            if ID not in ev_atom:\n",
    "                ev_atom.append(ID)\n",
    "\n",
    "            if line_index == len(ev_lines)-1:    # If at the last entry, finish clustering\n",
    "                \n",
    "                ev_no += 1\n",
    "                tf = ev_lines[line_index-1][0]\n",
    "\n",
    "                if tf == ti:\n",
    "                    dt = ev_lines[line_index-1][1]\n",
    "                else:\n",
    "                    dt = tf - ti\n",
    "\n",
    "                for atom in ev_atom:\n",
    "                    atom_flag = []    # List of flags for each atom\n",
    "                    \n",
    "                    for flag in ev_flag:\n",
    "                        ID = flag[2]\n",
    "                        \n",
    "                        if ID == atom:\n",
    "                            atom_flag.append(flag)\n",
    "                    \n",
    "                    for flag_index, flag in enumerate(atom_flag):\n",
    "                        \n",
    "                        ID = flag[2]\n",
    "                        Type = flag[3]\n",
    "                        ni = flag[4]\n",
    "                        nf = flag[5]\n",
    "                        CNi = flag[6]\n",
    "                        CNf = flag[7]\n",
    "                        dn = nf - ni\n",
    "                        \n",
    "                        if abs(dn) < 0.10*dn_avg:    # If normal does not change, print only if at the last flag\n",
    "                            \n",
    "                            if flag_index != len(atom_flag)-1:\n",
    "                                continue\n",
    "                            \n",
    "                            else:\n",
    "                                out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(ID, Type, ni, nf, CNi, CNf)\n",
    "                                out_cluster.write(out)\n",
    "                        \n",
    "                        elif abs(dn) > 0.90*dn_avg:    # If normal changes, print that flag and move on to the next atom\n",
    "                            out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(ID, Type, ni, nf, CNi, CNf)\n",
    "                            out_cluster.write(out)\n",
    "                            break\n",
    "\n",
    "                out = '%i\\t%i\\t%i\\n\\n'%(ev_no, tf, dt)    # Print the event number, time step, and duration\n",
    "                out_cluster.write(out)\n",
    "                    \n",
    "        # If a step is after 3 steps or more from the previously flagged step, it is part of a different event, so finish clustering\n",
    "        else:\n",
    "\n",
    "            ev_no += 1\n",
    "            tf = ev_lines[line_index-1][0]\n",
    "\n",
    "            if tf == ti:\n",
    "                dt = ev_lines[line_index-1][1]\n",
    "            else:\n",
    "                dt = tf - ti\n",
    "\n",
    "            for atom in ev_atom:\n",
    "                atom_flag = []\n",
    "                \n",
    "                for flag in ev_flag:\n",
    "                    ID = flag[2]\n",
    "                    \n",
    "                    if ID == atom:\n",
    "                        atom_flag.append(flag)\n",
    "                \n",
    "                for flag_index, flag in enumerate(atom_flag):\n",
    "                    \n",
    "                    ID = flag[2]\n",
    "                    Type = flag[3]\n",
    "                    ni = flag[4]\n",
    "                    nf = flag[5]\n",
    "                    CNi = flag[6]\n",
    "                    CNf = flag[7]\n",
    "                    dn = nf - ni\n",
    "                    \n",
    "                    if abs(dn) < 0.10*dn_avg:\n",
    "                        \n",
    "                        if flag_index != len(atom_flag)-1:\n",
    "                            continue\n",
    "                        \n",
    "                        else:\n",
    "                            out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(ID, Type, ni, nf, CNi, CNf)\n",
    "                            out_cluster.write(out)\n",
    "                    \n",
    "                    elif abs(dn) > 0.90*dn_avg:\n",
    "                        out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(ID, Type, ni, nf, CNi, CNf)\n",
    "                        out_cluster.write(out)\n",
    "                        break\n",
    "            \n",
    "            out = '%i\\t%i\\t%i\\n\\n'%(ev_no, tf, dt)\n",
    "            out_cluster.write(out)\n",
    "            \n",
    "            ev_atom = []    # Reset for the next event\n",
    "            ev_flag = []\n",
    "            \n",
    "            ev_flag.append(line)\n",
    "            ti = line[0] - line[1]\n",
    "            ID = line[2]\n",
    "            \n",
    "            if ID not in ev_atom:\n",
    "                ev_atom.append(ID)\n",
    "\n",
    "out_cluster.close()\n",
    "\n",
    "print('Event clustering done: %i ' %(ev_no) + 'event clusters identified > event_cluster.txt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Event localization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sitator.util import PBCCalculator\n",
    "pbc = PBCCalculator(lat)    # Periodic boundary condition\n",
    "\n",
    "cluster = open(file_cluster,'r')    # Load event clusters\n",
    "cluster_lines = cluster.readlines()\n",
    "cluster_lines = [line.split() for line in cluster_lines]\n",
    "\n",
    "file_ls = pwd + '/event_list.txt'    # List of localized events\n",
    "out_ls = open(file_ls,'w')\n",
    "\n",
    "out_ls.write(header1)\n",
    "out_ls.write(header2)\n",
    "\n",
    "for line_index, line in enumerate(cluster_lines):\n",
    "    if line and (line_index > 10):\n",
    "        if len(line) == 6:\n",
    "            line[0:2] = np.array(line[0:2]).astype(int)    # Atom ID, atom type\n",
    "            line[2:4] = np.array(line[2:4]).astype(float)    # Initial normal, final normal\n",
    "            line[4:6] = np.array(line[4:6]).astype(int)    # Initial coordination, final coordination\n",
    "        else:\n",
    "            line[0:3] = np.array(line[0:3]).astype(int)    # event number, time step, duration\n",
    "\n",
    "line_cont = 10    # Line number continuation\n",
    "N_splice = 0    # Number of events spliced\n",
    "while line_cont < len(cluster_lines)-2:\n",
    "    \n",
    "    ev_atom = []    # List of event atoms\n",
    "    for line_index, line in enumerate(cluster_lines):\n",
    "        if line and (line_index > line_cont):\n",
    "            if len(line) == 6:\n",
    "                ev_atom.append(int(line[0]))\n",
    "            else:\n",
    "                ev_no = line[0]\n",
    "                tf = line[1]\n",
    "                dt = line[2]\n",
    "                break\n",
    "    \n",
    "    ev_coord = []    # List of Cartesian coordinates for event atoms\n",
    "    for ID in ev_atom:\n",
    "        for atom in range(0,N_dump):\n",
    "            if xyz.finaldict[0][1]['id'][atom] == ID:\n",
    "                index = atom\n",
    "                x = xyz.finaldict[step[tf]][1]['xu'][index]\n",
    "                y = xyz.finaldict[step[tf]][1]['yu'][index]\n",
    "                z = xyz.finaldict[step[tf]][1]['z'][index]\n",
    "                r = np.array([x, y, z])\n",
    "                ev_coord.append(r)\n",
    "                break\n",
    "    ev_coord = np.array(ev_coord)\n",
    "        \n",
    "    ev_rij = pbc.pairwise_distances(ev_coord)    # Pairwise distance matrix for event atoms\n",
    "    ev_connectivity = (ev_rij <= 2.0*dr_avg)    # Connectivity matrix (T/F) within a distance cutoff\n",
    "    \n",
    "    from scipy.sparse.csgraph import connected_components\n",
    "    N_group, ls_group = connected_components(ev_connectivity)    # N_group = number of groups; ls_group = list of groups for each event atom\n",
    "    \n",
    "    if N_group == 1:    # No splicing needed\n",
    "        \n",
    "        for line_index, line in enumerate(cluster_lines):\n",
    "            if line and (line_index > line_cont):\n",
    "                if len(line) == 6:\n",
    "                    out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(line[0], line[1], line[2], line[3], line[4], line[5])\n",
    "                    out_ls.write(out)\n",
    "                else:\n",
    "                    out = '%i\\t%i\\t%i\\n\\n'%(ev_no, tf, dt)\n",
    "                    out_ls.write(out)\n",
    "                    line_cont = line_index\n",
    "                    break\n",
    "    \n",
    "    else:    # Splicing needed\n",
    "        \n",
    "        N_splice += 1\n",
    "        for group_index, group in enumerate(range(0,N_group)):\n",
    "            \n",
    "            ev_group = []    # List of event atoms for each group\n",
    "            for atom_index, atom_group in ls_group:\n",
    "                if atom_group == group:\n",
    "                    ev_group.append(ev_atom[atom_index])\n",
    "            \n",
    "            for line_index, line in enumerate(cluster_lines):    # Check whether group is z-active\n",
    "                if line and (line_index > line_cont):\n",
    "                    if len(line) == 6:\n",
    "                        if line[0] in ev_group:\n",
    "                            if abs(atom[2] - atom[3]) > 0.50*dn_avg:\n",
    "                                z_active = True\n",
    "                                break\n",
    "                    else:\n",
    "                        z_active = False\n",
    "                        break\n",
    "            \n",
    "            if z_active:    # Write only if group is z-active\n",
    "            \n",
    "                for line_index, line in enumerate(cluster_lines):\n",
    "                    if line and (line_index > line_cont):\n",
    "                        if len(line) == 6:\n",
    "                            if line[0] in ev_group:\n",
    "                                out = '\\t%i\\t%i\\t%f\\t%f\\t%i\\t%i\\n'%(line[0], line[1], line[2], line[3], line[4], line[5])\n",
    "                                out_ls.write(out)\n",
    "                        else:\n",
    "                            ev_no_new = str(ev_no) + '-' + str(group_index+1)\n",
    "                            out = '%s\\t%i\\t%i\\n\\n'%(ev_no_new, tf, dt)\n",
    "                            out_ls.write(out)\n",
    "                            if group_index == N_group-1:\n",
    "                                line_cont = line_index\n",
    "                            break                \n",
    "\n",
    "cluster.close()\n",
    "out_ls.close()\n",
    "\n",
    "print('Event localization done: %i ' %(N_splice) + 'events spliced for simultaneous events > event_list_local.txt\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Event classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_ls = pwd + '/event_list.txt'    # List of events\n",
    "ls = open(file_ls,'r')    # Load list of events\n",
    "ls_lines = ls.readlines()\n",
    "ls_lines = [line.split() for line in ls_lines]\n",
    "\n",
    "file_cls = pwd + '/event_class.txt'    # List of event classes\n",
    "out_cls = open(file_cls,'w')\n",
    "\n",
    "disclaimer = 'Intended only for estimation of the event statistics. Please confirm with visual inspection as needed.\\n'\n",
    "header = 'ID \\t Type \\t Site_i \\t Site_f \\n\\n'\n",
    "out_cls.write(disclaimer)\n",
    "out_cls.write(header)\n",
    "\n",
    "Lv1_1 = 'Isolated adatom'\n",
    "Lv1_2 = 'Edge/kink/vacancy'\n",
    "Lv1_3 = 'Deposit intralayer'\n",
    "Lv0_1 = 'Surface intralayer'\n",
    "Lv0_2 = 'Subdeposit'\n",
    "Lv_sub = 'Subsurface'\n",
    "\n",
    "for line_index, line in enumerate(ls_lines):\n",
    "    if line and (line_index > 9):\n",
    "        if len(line) == 6:\n",
    "            line[0:2] = np.array(line[0:2]).astype(int)    # Atom ID, atom type\n",
    "            line[2:4] = np.array(line[2:4]).astype(float)    # Initial normal, final normal\n",
    "            line[4:6] = np.array(line[4:6]).astype(int)    # Initial coordination, final coordination\n",
    "        else:\n",
    "            line[0:3] = np.array(line[0:3]).astype(int)   # Event #, time step, duration\n",
    "    \n",
    "line_cont = 9    # Line number continuation\n",
    "while line_cont < len(ls_lines)-2:\n",
    "    \n",
    "    ev_atom = []    # List of atom attributes for each event\n",
    "    for line_index, line in enumerate(ls_lines):\n",
    "        if line and (line_index > line_cont):\n",
    "            if len(line) == 6:\n",
    "                ev_atom.append(line)\n",
    "            else:\n",
    "                ev_no = line[0]\n",
    "                tf = line[1]\n",
    "                dt = line[2]\n",
    "                line_cont = line_index\n",
    "                break\n",
    "    \n",
    "    out = '%s%i\\n%s%i%s\\n%s%i%s\\n'%('Event #', ev_no, 'tf = ', tf, ' ps', 'dt = ', dt, ' ps')\n",
    "    out_cls.write(out)\n",
    "    \n",
    "    N = len(ev_atom)    # Number of active atoms\n",
    "    index_z = []    # Line indices of z-active atoms\n",
    "    index_aux = []    # Line indices of auxiliary atoms\n",
    "    \n",
    "    for atom_index, atom in enumerate(ev_atom):\n",
    "        if abs(atom[2] - atom[3]) > 0.50*dn_avg:\n",
    "            index_z.append(atom_index)\n",
    "        else:\n",
    "            index_aux.append(atom_index)\n",
    "    \n",
    "    N_z = len(index_z)    # Number of z-active atoms\n",
    "    N_aux = len(index_aux)    # Number of auxiliary atoms\n",
    "    \n",
    "    site = []    # Initial & final site & type matrix\n",
    "    for atom in ev_atom:\n",
    "        \n",
    "        ID = atom[0]\n",
    "        Type = atom[1]\n",
    "        n = [atom[2], atom[3]]\n",
    "        CN = [atom[4], atom[5]]\n",
    "        atom_site = [[],[]]\n",
    "        \n",
    "        for i in range(0,2):\n",
    "            \n",
    "            if n[i] < z_layer[N_slab-1]:\n",
    "                atom_site[i] = Lv_sub    # Subsurface\n",
    "            \n",
    "            elif n[i] == z_layer[N_slab-1]:\n",
    "                \n",
    "                if 11 <= CN[i] <= 12:\n",
    "                    atom_site[i] = Lv0_2    # Subdeposit\n",
    "                \n",
    "                elif 8 <= CN[i] <= 10:\n",
    "                    atom_site[i] = Lv0_1    # Surface intralayer\n",
    "            \n",
    "            elif n[i] > z_layer[N_slab-1]:\n",
    "                \n",
    "                if 8 <= CN[i] <= 12:\n",
    "                    atom_site[i] = Lv1_3    # Deposit intralayer\n",
    "                \n",
    "                elif 4 <= CN[i] <= 7:\n",
    "                    atom_site[i] = Lv1_2    # Edge/kink/vacancy\n",
    "            \n",
    "                elif 1 <= CN[i] <= 3:\n",
    "                    atom_site[i] = Lv1_1    # Isolated adatom\n",
    "        \n",
    "        site.append(atom_site)\n",
    "        out = '%i\\t%i\\t%s %s %s\\n'%(ID, Type, atom_site[0], '==>', atom_site[1])\n",
    "        out_cls.write(out)\n",
    "\n",
    "    # More than 5 active atoms = Indeterminate\n",
    "    if N > 5:\n",
    "        out = 'Visual inspection needed: More than 5 active atoms. Most likely conjoined events or intralayer shift.\\n'\n",
    "        out_cls.write(out)\n",
    "\n",
    "    # 1 z-active atom = Popout vs. insertion vs. step motion\n",
    "    elif N_z == 1:\n",
    "        \n",
    "        Type = ev_atom[index_z[0]][1]\n",
    "        \n",
    "        if Lv_sub in site[index_z[0]]:\n",
    "            out = 'Class = Interlayer transfer\\nAtom type = ' + str(Type) + '\\n'\n",
    "            out_cls.write(out)\n",
    "        \n",
    "        elif ((site[index_z[0]][0] == Lv0_1) or (site[index_z[0]][0] == Lv0_2)) and ((site[index_z[0]][1] == Lv1_1) or (site[index_z[0]][1] == Lv1_2) or (site[index_z[0]][1] == Lv1_3)):\n",
    "            out = 'Class = Popout\\nAtom type = ' + str(Type) + '\\n'\n",
    "            out_cls.write(out)\n",
    "        \n",
    "        elif ((site[index_z[0]][0] == Lv1_1) or (site[index_z[0]][0] == Lv1_2) or (site[index_z[0]][0] == Lv1_3)) and ((site[index_z[0]][1] == Lv0_1) or (site[index_z[0]][1] == Lv0_2)):\n",
    "            out = 'Class = Vacancy insertion\\nAtom type = ' + str(Type) + '\\n'\n",
    "            out_cls.write(out)\n",
    "        \n",
    "        elif ((site[index_z[0]][0] == Lv1_1) or (site[index_z[0]][0] == Lv1_2) or (site[index_z[0]][0] == Lv1_3)) and ((site[index_z[0]][1] == Lv1_1) or (site[index_z[0]][1] == Lv1_2) or (site[index_z[0]][1] == Lv1_3)):\n",
    "            \n",
    "            dn = ev_atom[index_z[0]][3] - ev_atom[index_z[0]][2]\n",
    "            \n",
    "            if N_aux != 0:\n",
    "                if dn > 0:\n",
    "                    out = 'Class = Exchange ascent\\nAtom type = ' + str(Type) + '\\n'\n",
    "                    out_cls.write(out)\n",
    "                else:\n",
    "                    out = 'Class = Exchange descent\\nAtom type = ' + str(Type) + '\\n'\n",
    "                    out_cls.write(out)\n",
    "            \n",
    "            else:\n",
    "                if dn > 0:\n",
    "                    out = 'Class = Hopping ascent\\nAtom type = ' + str(Type) + '\\n'\n",
    "                    out_cls.write(out)\n",
    "                else:\n",
    "                    out = 'Class = Hopping descent\\nAtom type = ' + str(Type) + '\\n'\n",
    "                    out_cls.write(out)\n",
    "    \n",
    "    # 2 z-active atoms = Popout vs. insertion vs. adatom exchange\n",
    "    elif N_z == 2:\n",
    "        \n",
    "        ID1 = ev_atom[index_z[0]][0]\n",
    "        ID2 = ev_atom[index_z[1]][0]\n",
    "        Type1 = ev_atom[index_z[0]][1]\n",
    "        Type2 = ev_atom[index_z[1]][1]\n",
    "        dn1 = ev_atom[index_z[0]][3] - ev_atom[index_z[0]][2]\n",
    "        dn2 = ev_atom[index_z[1]][3] - ev_atom[index_z[1]][2]\n",
    "\n",
    "        if Lv_sub in site[index_z[0]]:\n",
    "            out = 'Class = Interlayer transfer\\nAtom type = ' + str(Type1) + '\\n'\n",
    "            out_cls.write(out)\n",
    "            \n",
    "        elif Lv_sub in site[index_z[1]]:\n",
    "            out = 'Class = Interlayer transfer\\nAtom type = ' + str(Type2) + '\\n'\n",
    "            out_cls.write(out)\n",
    "        \n",
    "        elif (dn1 > 0) and (dn2 > 0):\n",
    "            \n",
    "            if ((site[index_z[0]][0] == Lv0_1) or (site[index_z[0]][0] == Lv0_2)):\n",
    "                out = 'Class = Popout\\nAtom type = ' + str(Type1) + '\\n'\n",
    "                out_cls.write(out)\n",
    "                \n",
    "            elif ((site[index_z[1]][0] == Lv0_1) or (site[index_z[1]][0] == Lv0_2)):\n",
    "                out = 'Class = Popout\\nAtom type = ' + str(Type2) + '\\n'\n",
    "                out_cls.write(out)\n",
    "            \n",
    "            else:\n",
    "                out = 'Visual inspection needed: Most likely simultaneous ascent.\\n'\n",
    "                out_cls.write(out)\n",
    "        \n",
    "        elif (dn1 < 0) and (dn2 < 0):\n",
    "            \n",
    "            if ((site[index_z[0]][1] == Lv0_1) or (site[index_z[0]][1] == Lv0_2)):\n",
    "                out = 'Class = Vacancy insertion\\nAtom type = ' + str(Type1) + '\\n'\n",
    "                out_cls.write(out)\n",
    "\n",
    "            elif ((site[index_z[1]][1] == Lv0_1) or (site[index_z[1]][1] == Lv0_2)):\n",
    "                out = 'Class = Vacancy insertion\\nAtom type = ' + str(Type2) + '\\n'\n",
    "                out_cls.write(out)\n",
    "            \n",
    "            else:\n",
    "                out = 'Visual inspection needed: Most likely simultaneous descent.\\n'\n",
    "                out_cls.write(out)\n",
    "        \n",
    "        else:\n",
    "            \n",
    "            if N_aux != 0:\n",
    "                out = 'Class = Indirect exchange\\nAtom type = ' + str(Type1) + '+' + str(Type2) + '\\n'\n",
    "                out_cls.write(out)\n",
    "            \n",
    "            else:    # Distance check\n",
    "                \n",
    "                for atom in range(0,N_dump1):\n",
    "                    if xyz1.finaldict[0][1]['id'][atom] == ID1:\n",
    "                        index1 = atom\n",
    "                    elif xyz1.finaldict[0][1]['id'][atom] == ID2:\n",
    "                        index2 = atom\n",
    "                \n",
    "                x1 = xyz1.finaldict[step[tf]][1]['xu'][index1]\n",
    "                y1 = xyz1.finaldict[step[tf]][1]['yu'][index1]\n",
    "                z1 = xyz1.finaldict[step[tf]][1]['z'][index1]\n",
    "                r1 = np.array([x1, y1, z1])\n",
    "                \n",
    "                x2 = xyz1.finaldict[step[tf]][1]['xu'][index2]\n",
    "                y2 = xyz1.finaldict[step[tf]][1]['yu'][index2]\n",
    "                z2 = xyz1.finaldict[step[tf]][1]['z'][index2]\n",
    "                r2 = np.array([x2, y2, z2])\n",
    "                \n",
    "                dr = np.subtract(r1,r2)\n",
    "                \n",
    "                if np.linalg.norm(dr) > 1.10*dr_avg:\n",
    "                    out = 'Class = Indirect exchange\\nAtom type = ' + str(Type1) + '+' + str(Type2) + '\\n' + 'Visual inspection advised: Missing auxiliary atom.\\n'\n",
    "                    out_cls.write(out)\n",
    "                else:\n",
    "                    out = 'Class = Direct exchange\\nAtom type = ' + str(Type1) + '+' + str(Type2) + '\\n'\n",
    "                    out_cls.write(out)\n",
    "    \n",
    "    # More than 2 z-active atoms = Interlayer exchange vs. indeterminate\n",
    "    elif N_z > 2:\n",
    "        \n",
    "        N_sub = 0\n",
    "        for index, i in enumerate(index_z):\n",
    "            \n",
    "            if Lv_sub in site[i]:\n",
    "                N_sub += 1\n",
    "                index_last = index\n",
    "            \n",
    "            if index == len(index_z)-1:\n",
    "                \n",
    "                if N_sub == 0:\n",
    "                    out = 'Visual inspection needed: More than two z-active atoms. Most likely conjoined events, or extra fluctuating atoms.\\n'\n",
    "                    out_cls.write(out)\n",
    "                \n",
    "                elif N_sub == 1:                    \n",
    "                    Type = ev_atom[index_last][1]\n",
    "                    out = 'Class = Interlayer transfer\\nAtom type = ' + str(Type) + '\\n'\n",
    "                    out_cls.write(out)\n",
    "\n",
    "                else:\n",
    "                    out = 'Class = Interlayer transfer\\nAtom type = Multi\\n'\n",
    "                    out_cls.write(out)\n",
    "\n",
    "    out_cls.write('\\n')\n",
    "    \n",
    "out_cls.close()\n",
    "ls.close()\n",
    "\n",
    "print('Event classification done > event_class.txt\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Event statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_cls = pwd + '/event_class.txt'    # List of event classes\n",
    "cls = open(file_cls,'r')    # Load list of event classes\n",
    "cls_lines = cls.readlines()\n",
    "\n",
    "file_ind = pwd + '/event_indeterminate.txt'    # Indeterminate events\n",
    "out_ind = open(file_ind,'w')\n",
    "\n",
    "file_stat = pwd + '/event_stat.txt'    # Event statistics\n",
    "out_stat = open(file_stat,'w')\n",
    "\n",
    "N_event = 0    # Number of events\n",
    "N_indeterminate = 0    # Number of indeterminate events\n",
    "for line_index, line in enumerate(cls_lines):\n",
    "    if 'Event' in line:\n",
    "        N_event += 1\n",
    "    elif 'Visual inspection needed:' in line:\n",
    "        N_indeterminate += 1\n",
    "        \n",
    "        end = line_index\n",
    "        for i in reversed(range(end)):\n",
    "            if 'Event' in cls_lines[i]:\n",
    "                start = i\n",
    "                break\n",
    "        \n",
    "        for i in range(start,end+1):\n",
    "            out_ind.write(cls_lines[i])\n",
    "        out_ind.write('\\n')\n",
    "        \n",
    "indeterminacy = 100*N_indeterminate/N_event\n",
    "header1 = 'Indeterminacy = %.2f' % (indeterminacy) + '% (' + str(N_indeterminate) + ' events)\\n'\n",
    "out_stat.write(header1)\n",
    "\n",
    "Class = ['Popout', 'Vacancy insertion', 'Exchange descent', 'Exchange ascent', 'Hopping descent', 'Hopping ascent', 'Indirect exchange', 'Direct exchange', 'Interlayer transfer']\n",
    "if N_type == 1:\n",
    "    hist = np.zeros((9,), dtype=int)\n",
    "else:\n",
    "    hist = np.zeros((21,), dtype=int)\n",
    "\n",
    "for line_index, line in enumerate(cls_lines):    # Count number of events in each class\n",
    "    if 'Class' in line:\n",
    "        \n",
    "        for C_index, C in enumerate(Class):\n",
    "            if C in line:\n",
    "                \n",
    "                for Type_index, Type in enumerate(range(1,N_type+1)):\n",
    "                    \n",
    "                    if C_index <= 5:\n",
    "                        if str(Type) in cls_lines[line_index+1]:\n",
    "                            hist[2*C_index + Type_index] += 1\n",
    "                            \n",
    "                    elif 6 <= C_index <= 7:\n",
    "                        \n",
    "                        Atom_type = str(Type) + '+' + str(Type)\n",
    "                        if Atom_type in cls_lines[line_index+1]:\n",
    "                            hist[2*C_index + Type_index + (C_index-6)] += 1\n",
    "                        \n",
    "                        if Type_index == 1:\n",
    "                            Atom_type1 = str(Type-1) + '+' + str(Type)\n",
    "                            Atom_type2 = str(Type) + '+' + str(Type-1)\n",
    "                            if (Atom_type1 in cls_lines[line_index+1]) or (Atom_type2 in cls_lines[line_index+1]):\n",
    "                                hist[2*C_index + Type_index + (C_index-6) + 1] += 1\n",
    "                    \n",
    "                    else:\n",
    "                        \n",
    "                        if str(Type) in cls_lines[line_index+1]:\n",
    "                            hist[2*C_index + Type_index + 2] += 1\n",
    "                        \n",
    "                        if Type_index == 1:\n",
    "                            if 'Multi' in cls_lines[line_index+1]:\n",
    "                                hist[-1] += 1\n",
    "                    \n",
    "total = np.sum(hist)\n",
    "header2 = str(total) + ' out of ' + str(N_event) + ' events classified.\\n\\n'\n",
    "out_stat.write(header2)\n",
    "\n",
    "header3 = ['Event class', 'Atom type', '# of occurrence', '% frequency']\n",
    "types = ['%%%is', '%%%is', '%%%ii', '%%%i.2f']\n",
    "colwidth = max(max(len(h) for h in header3), max(len(c) for c in Class))\n",
    "out_stat.write((('\\t'.join([\"%%%is\" % colwidth] * len(header3))) % tuple(header3)) + \"\\n\")\n",
    "\n",
    "for C_index, C in enumerate(Class):\n",
    "    for Type_index, Type in enumerate(range(1,N_type+1)):\n",
    "        \n",
    "        if C_index <= 5:\n",
    "            number = hist[2*C_index + Type_index]\n",
    "            freq = float(100*number/total)\n",
    "            out = ('\\t'.join(t % colwidth for t in types)) % (C, str(Type), number, freq)\n",
    "            out += '\\n'\n",
    "            out_stat.write(out)\n",
    "            \n",
    "        elif 6 <= C_index <= 7:\n",
    "            \n",
    "            Atom_type = str(Type) + '+' + str(Type)\n",
    "            number = hist[2*C_index + Type_index + (C_index-6)]\n",
    "            freq = float(100*number/total)\n",
    "            out = ('\\t'.join(t % colwidth for t in types)) % (C, Atom_type, number, freq)\n",
    "            out += '\\n'\n",
    "            out_stat.write(out)\n",
    "\n",
    "            if Type_index == 1:\n",
    "                Atom_type = str(Type-1) + '+' + str(Type)\n",
    "                number = hist[2*C_index + Type_index + (C_index-6) + 1]\n",
    "                freq = float(100*number/total)\n",
    "                out = ('\\t'.join(t % colwidth for t in types)) % (C, Atom_type, number, freq)\n",
    "                out += '\\n'\n",
    "                out_stat.write(out)\n",
    "        \n",
    "        else:\n",
    "            \n",
    "            number = hist[2*C_index + Type_index + 2]\n",
    "            freq = float(100*number/total)\n",
    "            out = ('\\t'.join(t % colwidth for t in types)) % (C, str(Type), number, freq)\n",
    "            out += '\\n'\n",
    "            out_stat.write(out)\n",
    "\n",
    "            if Type_index == 1:\n",
    "                Atom_type = 'Multi'\n",
    "                number = hist[-1]\n",
    "                freq = float(100*number/total)\n",
    "                out = ('\\t'.join(t % colwidth for t in types)) % (C, Atom_type, number, freq)\n",
    "                out += '\\n'\n",
    "                out_stat.write(out)\n",
    "\n",
    "out_stat.close()\n",
    "out_ind.close()\n",
    "cls.close()\n",
    "\n",
    "print('Event statistics done > event_stat.txt\\n')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
